Model Sensitivity

Estimating The Forward Electricity Curve In Brazil With A Model Of Two Agents Using Contracts By Difference And ECP_G Function

Authors: Felipe Van de Sande Araujo, Cristina Spineti Luz, Leonardo Lima Gomes, Luís Eduardo Teixeira Brandão

Abstract: The development of simple and effective mechanisms to estimate the value of the forward curve of power could enable market participants to better price hedging or speculative positions. This could in turn provide transparency in future price definition to all market participants and lead to more safety and liquidity in the market for electricity futures and power derivatives. This work presents a model for two market participants, a buyer and a seller of a contract for difference on the future spot price of electricity in southwest Brazil. It is shown that this model is representative of all market participants that have exposure to the future price of power. Each participant’s utility function is modeled using a Generalized Extended CVaR Preference (ECP_G) and the market equilibrium is obtained through the minimization of the quadratic difference between the certainty equivalent of both agents. The results are compared with prediction of the future spot price of power made by market specialists and found to yield reasonable results when using out of sample data.

This work presents the calculations done for the article referenced above. The following calculations show the sensitivity analysis for the optimized parameters which were obtained in the optimization step (link) and the risk-free return rate used in the validation step (link). All the code was run in RStudio using the version of the software below.

Software version

R.version

               _                           
platform       x86_64-w64-mingw32          
arch           x86_64                      
os             mingw32                     
system         x86_64, mingw32             
status                                     
major          4                           
minor          0.1                         
year           2020                        
month          06                          
day            06                          
svn rev        78648                       
language       R                           
version.string R version 4.0.1 (2020-06-06)
nickname       See Things Now              

Setting of the environment

# Set plot font
windowsFonts(A = windowsFont("Times New Roman")) 

Local functions definition

The first function is the determination of the ECPG phi which is the price which, when added to the contract price, provides the certainty equivalent for the agents. This is clarified in the text.

ECPG.phi <- function(parameters, dataset, auxList, buyer=TRUE) {
  if (any(parameters < 0) || any(parameters > 1)) return(NA)
  sumLambda = parameters['Lambda1'] + parameters['Lambda2']
  if (sumLambda > 1) return(NA)
  Lambda0 = 1 - sumLambda
  if (Lambda0 < 0) return(NA)
  
  if (buyer) {
    mean.pld = auxList$meanPLD
    VaR0 = auxList$VaR0
    VaR1 = auxList$VaR1
  } else {
    dataset = -dataset
    mean.pld = -auxList$meanPLD
    VaR0 = auxList$VaR1
    VaR1 = auxList$VaR0
  }
  
  VaRAlpha1 = apply(dataset, 1, quantile, (1-parameters['Alpha1']))
  VaRAlpha2 = apply(dataset, 1, quantile, (1-parameters['Alpha2']))
  
  CVaR1 = CVaR2 = NA
  
  CVaR1.data = dataset
  CVaR1.data[dataset > VaRAlpha1] <- NA
  
  CVaR2.data = dataset
  CVaR2.data[dataset > VaRAlpha2] <- NA
  
  CVaR1 = apply(CVaR1.data, 1, mean, na.rm=TRUE)
  CVaR2 = apply(CVaR2.data, 1, mean, na.rm=TRUE)
  
  Sum1 = parameters['Lambda1']*VaRAlpha1 + parameters['Lambda2']*VaRAlpha2
  Sum2 = parameters['Lambda1']*(CVaR1 - VaRAlpha1) + parameters['Lambda2']*(CVaR2 - VaRAlpha2)
  
  Limit1 = Lambda0*VaR0 + Sum1
  Limit2 = Lambda0*VaRAlpha1 + Sum1
  Limit3 = Lambda0*VaRAlpha2 + Sum1
  Limit4 = Lambda0*VaR1 + Sum1
  
  ECPG = Lambda0*mean.pld + parameters['Lambda1']*CVaR1 + parameters['Lambda2']*CVaR2
  
  nIter = nrow(dataset)
  Lambda1 = rep(parameters['Lambda1'], nIter)
  Lambda2 = rep(parameters['Lambda2'], nIter)
  
  if (any(ECPG>=Limit2)) Lambda1[ECPG>=Limit2] = 0
  if (any(ECPG>=Limit3)) Lambda2[ECPG>=Limit3] = 0
  
  Q = Lambda0 + Lambda1/(1-parameters['Alpha1']) + Lambda2/(1-parameters['Alpha2'])
  
  A = Sum2 + Lambda0*mean.pld + 
    Lambda1*VaRAlpha1/(1-parameters['Alpha1']) +
    Lambda2*VaRAlpha2/(1-parameters['Alpha2'])
  
  Eq = A/Q
  
  if (any(is.na(Eq))) return(NA)
  return(Eq)
}

We will use a list with all the parameters and data which will be called ECPG_Object, in order to facilitate the next function calls. The next function will create the list.

createECPGObject <- function(dataset, discountVector, forwardPrice, buyerParam, sellerParam) {
  VaR0 = apply(dataset, 1, quantile, 1)
  VaR1 = apply(dataset, 1, quantile, 0)
  meanPLD = apply(dataset, 1, mean)
  
  ECPGObject = list(
    dataset=dataset,
    discountVector=discountVector,
    forwardPrice=forwardPrice,
    buyerParam=buyerParam,
    sellerParam=sellerParam,
    VaR0=VaR0,
    VaR1=VaR1,
    meanPLD=meanPLD
  )
  return(ECPGObject)
}

Then we will define a function to obtain the equilibrium price using the calculated phi for each agent. The distinction between the agents must be made because the signal of the price data changes between them, because of the Contract for Differences revenue equations.

ECPG.equilibrium <- function (ECPGobject) {
  auxList = list(
    VaR0 = ECPGobject$VaR0,
    VaR1 = ECPGobject$VaR1,
    meanPLD = ECPGobject$meanPLD
  )
  
  buyer.phi = ECPG.phi(ECPGobject$buyerParam, ECPGobject$dataset, auxList)
  seller.phi = -ECPG.phi(ECPGobject$sellerParam, ECPGobject$dataset, auxList, buyer=FALSE)
  
  equilibrium.prices = (buyer.phi + seller.phi)/2
  
  adjusted.prices = equilibrium.prices/ECPGobject$discount
  
  averagePrice = mean(adjusted.prices)
  
  ECPGobject$buyerPhi = buyer.phi
  ECPGobject$sellerPhi = seller.phi
  ECPGobject$equilibriumPrices = equilibrium.prices
  ECPGobject$adjustedPrices=adjusted.prices
  ECPGobject$averagePrice=averagePrice
  
  return(ECPGobject)
}

Next function is used in a recursive call to calculate sensitivity.

sensitivity.foo <- function(i, arg, output, combination.matrix, ECPGobject) {
  
  ECPGobject[[arg]][TRUE] = combination.matrix[i, ]
  
  return(ECPG.equilibrium(ECPGobject)[[output]])
}

And finally a plotting function for the ECPG_Object.

sensitivity.prettyPlot <- function(plot.matrix, xlabel="", ylim=c(0,250), lwd = 2, legendTitle, legendPos="bottomright", inset=0.05, ncol=1, cex=1, legCex=0.9) {
  Nseries = ncol(plot.matrix)
  viridis.pallette = viridis(Nseries)
  line.types = rep(1, Nseries)
  
  Nrows = 1:nrow(plot.matrix)
  labels=row.names(plot.matrix)
  
  par(family="A", cex=cex)
  matplot(Nrows, plot.matrix, type='l', xlab=xlabel, ylab='Energy Price (R$/MWh)', ylim=ylim, lty=line.types, col=viridis.pallette, lwd=lwd, axes=F)
  axis(2)
  axis(side=1,at=Nrows,labels=labels)
  grid (NULL,NULL, lty = 6, col = "grey")
  legend(x=legendPos, inset=inset, legend=colnames(plot.matrix), horiz=FALSE, lty=line.types, col=viridis.pallette, lwd=lwd, cex=legCex, ncol=ncol, title=legendTitle)
}

Local variables definition

These variables will be used for all sensitivity tests. The optimization parameters are obtained from the training step.

 
buyer.param = c(Lambda1 = 0.5712411, Lambda2 = 0.1400453, Alpha1 = 0.5845306, Alpha2 = 0.5995087)
seller.param = c(Lambda1 = 0.3483552, Lambda2 = 0.5599081, Alpha1 = 0.7256726, Alpha2 = 0.9737610)

Variables for January A+1 analysis

Create a vector with discount rates to adjust the equilibrium prices with the risk-free return rate. Risk-free return rate is obtained as an approximation of the SELIC national rate.

# Set risk free rates
anualRate = 0.05
monthlyRate = (1+anualRate)^(1/12)-1
temp.discount.vector <- rep(NA, 23)
for (i in 1:23){
  temp.discount.vector[i] <- (1+monthlyRate)^(i-1)
}
januaryA1.discountVector = temp.discount.vector[12:23]

Input the forward price for electricity obtained from DCide Energia.

januaryA1.forwardPrice = 199.32

Loading the January A+1 data

Loading future spot price (PLD) data.

januaryA1.original.data = read.csv(file='./Data/DataJanuaryA1.csv', col.names=paste(month.abb, 2020, sep="-"))
januaryA1.data = t(januaryA1.original.data)

Create a list with all variables and data to pass on to functions.

januaryA1.ECPG_object = createECPGObject(dataset = januaryA1.data, buyerParam = buyer.param, sellerParam = seller.param, discountVector = januaryA1.discountVector, forwardPrice = januaryA1.forwardPrice)

Get the original results for comparison.

januaryA1.Results = ECPG.equilibrium(januaryA1.ECPG_object)

Sensitivity for Lambda 1 and Lambda 2

The first analysis is regarding \(\lambda_1\) and \(\lambda_2\) parameters for each agent. Those are the auxiliary variables set for this task.

lambdas.vec = seq(0, 0.45, by=0.05)
lambdas.matrix = as.matrix(expand.grid(lambdas.vec, lambdas.vec))
colnames(lambdas.matrix) = c('Lambda1', 'Lambda2')
nLambdas = nrow(lambdas.matrix)

From the range set above a combination matrix is created with the buyer parameters.

buyer.lambdas.matrix = cbind(lambdas.matrix, Alpha1=rep(buyer.param['Alpha1'], nLambdas), Alpha2=rep(buyer.param['Alpha2'], nLambdas))

Next, call the sensitivity function and create a matrix of the results for the buyer.

lambdaBuyer.sensitivity = sapply(1:nLambdas,
                                 sensitivity.foo,
                                 arg = "buyerParam",
                                 output = "averagePrice",
                                 combination.matrix = buyer.lambdas.matrix,
                                 januaryA1.ECPG_object)

lambdaBuyer.sensiMatrix = matrix(lambdaBuyer.sensitivity,
                                 ncol=length(lambdas.vec),
                                 byrow = TRUE,
                                 dimnames = list(lambdas.vec, lambdas.vec))

Do the same for the seller.

seller.lambdas.matrix = cbind(lambdas.matrix, Alpha1=rep(seller.param['Alpha1'], nLambdas), Alpha2=rep(seller.param['Alpha2'], nLambdas))

lambdaSeller.sensitivity = sapply(1:nLambdas,
                                  sensitivity.foo,
                                  arg = "sellerParam",
                                  output = "averagePrice",
                                 combination.matrix = seller.lambdas.matrix,
                                 januaryA1.ECPG_object)

lambdaSeller.sensiMatrix = matrix(lambdaSeller.sensitivity,
                                 ncol=length(lambdas.vec),
                                 byrow = TRUE,
                                 dimnames = list(lambdas.vec, lambdas.vec))

The result are plotted in an interactive widget using plotly.

require(plotly)

Loading required package: plotly
Loading required package: ggplot2
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: 㤼㸱plotly㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

ax.xy <- list(ticketmode = 'array', ticktext = as.character(lambdas.vec), tickvals=0:9)
axz <- list(title = "Power Price (R$/MWh)")

Analysis of \(\lambda_1\) and \(\lambda_2\) for the buyer.

fig.Buyer <- plot_ly(width=700, height=500) %>%
  add_surface(z = lambdaBuyer.sensiMatrix) %>%
  layout(scene = list(xaxis=c(ax.xy, title = "Buyer Lambda1"),yaxis=c(ax.xy, title = "Buyer Lambda2"),zaxis=axz))

fig.Buyer

Analysis of \(\lambda_1\) and \(\lambda_2\) for the seller.

fig.Seller <- plot_ly(width=700, height=500) %>%
  add_surface(z = lambdaSeller.sensiMatrix) %>%
  layout(scene = list(xaxis=c(ax.xy, title = "Seller Lambda1"),yaxis=c(ax.xy, title = "Seller Lambda2"),zaxis=axz))

fig.Seller

Sensitivity for Lambda 1 and Alpha 1 in ECP

The following test is regarding both \(\alpha\) and \(\lambda\) parameters for each agent using ECP measure instead of ECP_G. Those are the auxiliary variables set for this task.

lambdaAlpha.vec = seq(0, 0.9, by=0.1)
lambdaAlpha.matrix = as.matrix(expand.grid(lambdaAlpha.vec, lambdaAlpha.vec))
colnames(lambdaAlpha.matrix) = c('Lambda1', 'Alpha1')
nLambdaAlpha = nrow(lambdaAlpha.matrix)

From the range set above a combination matrix is created with the buyer parameters.

buyer.lambdaAlpha.matrix = cbind(lambdaAlpha.matrix, Lambda2=rep(0, nLambdaAlpha), Alpha2=rep(buyer.param['Alpha2'], nLambdaAlpha))
buyer.lambdaAlpha.matrix = buyer.lambdaAlpha.matrix[, c('Lambda1', 'Lambda2', 'Alpha1', 'Alpha2')]

Next, call the sensitivity function and create a matrix of the results for the buyer.

lambdaAlphaBuyer.sensitivity = sapply(1:nLambdaAlpha,
                                 sensitivity.foo,
                                 arg="buyerParam",
                                 output = "averagePrice",
                                 combination.matrix = buyer.lambdaAlpha.matrix,
                                 within(januaryA1.ECPG_object, buyerParam[TRUE] <- rep(0,4)))

lambdaAlphaBuyer.sensiMatrix = matrix(lambdaAlphaBuyer.sensitivity,
                                 ncol=length(lambdaAlpha.vec),
                                 byrow = TRUE,
                                 dimnames = list(lambdaAlpha.vec, lambdaAlpha.vec))

And the same is done for the seller.

seller.lambdaAlpha.matrix = cbind(lambdaAlpha.matrix, Lambda2=rep(0, nLambdaAlpha), Alpha2=rep(seller.param['Alpha2'], nLambdaAlpha))
seller.lambdaAlpha.matrix = seller.lambdaAlpha.matrix[, c('Lambda1', 'Lambda2', 'Alpha1', 'Alpha2')]

lambdaAlphaSeller.sensitivity = sapply(1:nLambdaAlpha,
                                  sensitivity.foo,arg="sellerParam",
                                  output = "averagePrice",
                                  combination.matrix = seller.lambdaAlpha.matrix,
                                  within(januaryA1.ECPG_object, sellerParam[TRUE] <- rep(0,4)))

lambdaAlphaSeller.sensiMatrix = matrix(lambdaAlphaSeller.sensitivity,
                                  ncol=length(lambdaAlpha.vec),
                                  byrow = TRUE,
                                  dimnames = list(lambdaAlpha.vec, lambdaAlpha.vec))

The results are visualized with our own plot function. Library viridis is called to provide color pallette for colorblind visualization.

require(viridis)

Loading required package: viridis
Loading required package: viridisLite

Alpha-Lambda sensitivity for the buyer.

sensitivity.prettyPlot(lambdaAlphaBuyer.sensiMatrix, xlabel=expression(paste("Buyer ", alpha, 1)), ylim=c(170, 205), ncol=2, cex=1.2, legCex=0.7, lwd=2, legendTitle=expression(paste("Buyer ", lambda, 1)))

Alpha-Lambda sensitivity for the seller.

sensitivity.prettyPlot(lambdaAlphaSeller.sensiMatrix, xlabel=expression(paste("Seller ", alpha, 1)), ylim=c(0, 180), legendTitle=expression(paste("Seller ", lambda, 1)), legendPos = "topleft", ncol=5, inset=0.02, cex=1.2, lwd=2, legCex=0.7)

Sensitivity for the risk-free rate of return in January A+2 results

This section tests the result for both the average price and the adjusted prices obtained by the model with different levels of the risk-free rate.

Define Variables for January A+2 validation

Create a vector with discount rates to adjust the equilibrium prices with the risk-free return rate. At this point all values are set to 1.

januaryA2.discountVector = rep(1, 12)

Input the forward price for electricity obtained from DCide Energia.

januaryA2.forwardPrice = 179.26

Loading the January A+2 data

Loading future spot price (PLD) data.

januaryA2.original.data = read.csv(file='./Data/DataJanuaryA2.csv', col.names=paste(month.abb, 2021, sep="-"))
januaryA2.data = t(januaryA2.original.data)

Create a list with all variables and data to pass on to functions.

januaryA2.ECPG_object = createECPGObject(dataset = januaryA2.data, buyerParam = buyer.param, sellerParam = seller.param, discountVector = januaryA2.discountVector, forwardPrice = januaryA2.forwardPrice)

Prepare variables for the sensitivity analysis

Set auxiliary variables.

annualRates.vec = seq(0, 0.09, by=0.01)
monthlyRates.vec = (1+annualRates.vec)^(1/12)-1
nRates = length(monthlyRates.vec)

Create a discount matrix to be passed to the sensitivity function.

discount.matrix <- matrix(NA, ncol=35, nrow=nRates)
for (ii in 1:nRates) {
  for (i in 1:35){
    discount.matrix[ii, i] <- (1+monthlyRates.vec[ii])^(i-1)
  }
}
discount.matrix = discount.matrix[, 24:35]
rownames(discount.matrix) = annualRates.vec

Apply the sensitivity function to obtain the average price for each risk-free rate.

rates.sensitivity = sapply(1:nRates,
                            sensitivity.foo,
                            arg="discountVector",
                           output = "averagePrice", 
                            combination.matrix = discount.matrix,
                            januaryA2.ECPG_object)
names(rates.sensitivity) = annualRates.vec

Plot the results

plot(annualRates.vec, rates.sensitivity, ylim=c(0,200), ylab="Energy Price (R$/MWh)", xlab="Annual Risk-Free Rate", col="steelblue")

Apply the sensitivity function to obtain the adjusted prices for each risk-free rate and arrange the results in a matrix.

rates.adjustedPrices = sapply(1:nRates,
                           sensitivity.foo,
                           arg="discountVector",
                           output = "adjustedPrices", 
                           combination.matrix = discount.matrix,
                           januaryA2.ECPG_object)

rates.sensitivity.matrix = matrix(rates.adjustedPrices,
                                  ncol = nRates,
                                  byrow = FALSE,
                                  dimnames = list(rownames(januaryA2.ECPG_object$dataset), annualRates.vec))

Plot the results

sensitivity.prettyPlot(rates.sensitivity.matrix, ylim=c(0, 205), ncol=5, cex=1.2, legCex=0.7, lwd=2, legendTitle="Risk-Free Rate aa.")

Check if results are internally consistent.

januaryA2.Results = ECPG.equilibrium(within(januaryA2.ECPG_object, discountVector <- discount.matrix["0.05", ]))
januaryA2.Results$averagePrice == rates.sensitivity["0.05"]

0.05 
TRUE 

See also:

Model Training (link)

Model Validation (link)